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1 Introduction 

1.1 Ferromagnetism: whither theory? 

Metallic ferromagnetism has since long been a controversial subject in solid-state theory. 
From early on, two schools of thought have been competing. The first, based on ideas 
of Slater, Stoner, and Wohlfarth, turns to the band aspects of metallic ferromagnets 
such as iron, cobalt, and nickel PP. The itineracy of the conduction electrons, including 
those of the open 3d shell, is considered to dominate the electron-electron interaction. 
The second school, going back to van-Vleck |2j and Gutzwiller j^, stresses the localized 
character of the 3c? electrons in transition metals and their compounds; local magnetic 
moments are indeed manifest in magnetic insulators which are often oxides and halides 
of transition metals. In the van-Vleck-Gutzwiller school, local moments are seen as an 
important source of magnetism, common to both ferromagnetic transition metals and 
their insulating compounds. Therefore, the van-Vleck-Gutzwiller school emphasizes the 
importance of correlations as a consequence of the electron-electron interactions. 



As the Slater-Stoner-Wohlfarth school considers electronic correlations to be of minor 
importance, it treats the electron-electron interaction in simple approximations such as 
Hartree-Fock theory, see p. Later, the density- functional theory (DFT) in its spin- 
dependent versions (SDFT) took over the role of Hartree-Fock or Xa schemes of the Slater 
school. For more than two decades, SDFT energy-band calculations by various groups 
have provided reasonably accurate results for quite a number of properties of metallic 
ferromagnets, e.g., magnetic moments and the overall shape of the Fermi surface. An 
early problem was the observation that SDFT in standard local-density approximation 
(LDA) does not give the bcc lattice as the most stable crystal structure for iron, instead 
fee Fe is found to be stable. This problem was overcome by the introduction of generalized 
gradient expansions (GGA) of the density pi. 

Until recently the van-Vleck-Gutzwiller school could not provide much more than 
intuitive arguments for its case because ferromagnetism is seen as a genuine, difficult 
many-body problem. For decades, only simplistic models could be treated with some 
reliability, e.g., the one-band Hubbard model, which, however, exhibits ferromagnetism 
only for an extremely large Hubbard interaction or a pathological density of state 
In the past, despite some partial successes jH], correlated-electron theory could not de- 
scribe real materials in much detail. Over the last decade new theoretical tools have been 
developed which try to combine the merits of density-functional theory with those of 
advanced many-body approaches, e.g., the LDA-GW approximation |7j. A general prob- 
lem of correlated-electron theories is the construction of suitable model Hamiltonians. 
Progress in this direction has been made recently, e.g., the 'down-folding' scheme permits 
a systematic construction of Hamiltonians in a reduced Hilbert space [Hj. The one-particle 
dynamics of models for correlated lattice electrons can be determined, in principle, within 
the dynamical mean-field theory (DMFT) ^ ^1 which becomes exact in the limit 
of an infinite number of nearest neighbors (coordination number Z — > oo). The 'static 
approximation' to DMFT is the 'LDA+[/' scheme which has been applied to a variety of 
correlated-electron systems [T^ ITHl IT^ . Moreover, DMFT has been extended in various 
ways, e.g., a combination of the LDA-GW approximation and the DMFT has been put 
forward jT^I. Therefore, correlated-electron theories seem to approach a level of quality 
where a faithful comparison with experiments becomes possible. 

In our group, we have extended the original Gutzwiller variational scheme to the 
full multi-band problem. Our theory provides total energies and magnetic moments of 
competing magnetic and non-magnetic phases and the single-particle excitation energies 
within Fermi-liquid theory. Thereby, our theory describes the dispersion of the conduction 
and valence bands and the Fermi surface as accessible, e.g., by angle-resolved photo- 
emission spectroscopy (ARPES), and the predictions of our Gutzwiller correlated-electron 
theory can be tested against experiments on transition metals and its compounds. 

1.2 The nickel problem 

In this article we address ferromagnetic nickel as the prototype of band ferromagnetic 
materials. On the experimental side, one of the reasons why nickel has been in the focus 
of interest, much more than iron or cobalt, is the relative ease to grow large single crystals. 



Consequently, a lot of experimental investigations have been carried out on nickel in order 
to determine its magnetic properties and electronic structure. For a review, we refer to 
the articles [H CHI 113 HE] • 

On the theoretical side, of all the iron group magnetic metals nickel is the most cele- 
brated case of discrepancies between SDFT predictions and experimental results. Early 
on, de-Haas-van-Alphen data indicated the presence of only one hole ellipsoid in the 
Fermi surface of the minority spin bands, located around the X-point of the Brillouin 
zone ^ni- In contrast, energy band calculations using SDFT predict two hole ellipsoids 
around the X-point, the Xs-state of pure d(t2g) orbital character, and the X2-state of 
pure d{eg) orbital character. Concomitantly, neutron scattering experiments revealed 
that the orbital character of the nickel magnetic moment exhibits less eg admixture than 
predicted by energy band calculations. 

ARPES measurements confirmed the de-Haas-van-Alphen results: the minority X2- 
state of pure d{eg) orbital character was found to lie slightly below the Fermi level. Fur- 
thermore, the photoemission data differ considerably from SDFT over the whole range 
of occupied 3d bands. Most importantly, the width of the occupied part of the 3d 
bands is approximately W = 3.3 eV [211 122], whereas SDFT gives VTsdft = 4.5 eV, or 
larger ^[ . ARPES studies also reveal big discrepancies between SDFT and experiment in 
the exchange splittings of majority and minority spin bands near the Fermi energy. The 
SDFT results give a rather isotropic exchange splitting of about 650 meV to 750 meV. 
whereas photoemission data show small and highly anisotropic exchange splittings be- 
tween 160 meV for pure 3c?(eg)-states such as X2 and 330 meV for pure 3(i(t2g)-states such 
as S3 near the X-point. 

Another discrepancy concerns the position of the L2'-state. Experimentally, it is found 
approximately 1 eV below the Fermi energy E-p whereas the SDFT puts it to about 0.3 eV 
below Ep. This state is special as it represents the most binding, pure 4p-state 
An incorrect 4p-level implies that DFT slightly underestimates the partial density of the 
4p electrons and, thus, the 3d-hole density. 

Lastly, the size and the direction of the magnetic moment pose a big problem to SDFT. 
Early experiments using ferromagnetic electron resonance and the Einstein-de-Haas effect 
have produced the values g = 2.183 and g' = 1.835 which result in an orbital moment 
of /Xorb = 0.0507/xb per atom for nickel jT]. In general, the values of the orbital part 
of the magnetic moments in band ferromagnets are not very well reproduced in SDFT. 
Iron and cobalt show large deviations [23; even in relativistic SDFT calculations pB] . 
Moreover, the magneto-crystalline anisotropy energies are not correct, i.e., they come 
with the wrong sign in SDFT. As shown by Aubert and Gersdorf [2ZII2H]! the value of the 
dominant anisotropy constant is Ki = — 8.8/ieV, i.e., the (111) direction is the easy axis 
in nickel. For a (001) magnetic moment, the energy is higher with an energy difference 
£^aniso = Eui — Eqqi ^ —3 /ieV per atom. In contrast, it is -Eaniso > in SDFT calculations, 
e.g., by Daalderoop j2S] so that Ki > 0. The SDFT results for cobalt and iron are equally 
disappointing. At low temperatures, the complex behavior of the magnetic anisotropy 
as a function of temperature and the magnetic field was interpreted by Gersdorf j2Z] to 
be caused by a shift of the electronic state (001) X21 to about 3meV above the Fermi 
energy. The famous second hole ellipsoid around the X21 indeed exists |27j as long as the 



magnetic moment has an angle of less than 18 degrees with the (001) direction. However, 
the states (100) X21 and (010) X21 remain below E-p, about 33meV j2Z|. 

All these discrepancies could not be overcome by SDFT modifications such as the 
GGA variants Apparently, a substantial improvement of SDFT results for nickel 
along these lines is presently not within reach. Recent correlated-electron theories remove 
the SDFT shortcomings only partially. For example, with LDA+?7-type calculations it 
seems to be possible to adjust to experiment the orbital moment, the magnetic anisotropy, 
and the Fermi surface topology [131111; however, other quantities, especially the entire 
quasi-particle band structure, remain elusive within these correlated-electron theories. 

1.3 Scope of this work 

In this article we show that our Gutzwiller theory resolves all of the discrepancies described 
in Section 11.21 Gutzwiller theory reproduces the experimental Fermi surface topology, 
the quasi-particle band structure with the correct width and the anisotropic exchange 
splittings, and the magnetic anisotropy. This success for nickel corroborates the early 
views by van-Vleck and Gutzwiller: the inclusion of the electron-electron interaction on 
a many-body level is essential for a proper understanding of itinerant ferromagnetism. 

We structure our article as follows. In Section |21 we introduce our general Gutzwiller- 
correlated wave functions for multi-band Hubbard models. In Section El we present ana- 
lytical results in the limit of infinite coordination number which we apply to nickel. The 
remaining Sections EHZl are devoted to an improved, concise derivation of our theory. A 
short outlook closes our presentation. 

The class of variational wave functions for which our formulae apply goes much be- 
yond our previous investigations [^lllSniE]- First, the Gutzwiller correlator now employs 
projections onto eigenstates of an effective atomic Hamiltonian which makes the varia- 
tional Hilbert space much more fiexible with respect to local correlations. Second, our 
variational states cover the case of superconductivity in correlated multi-band lattice sys- 
tems. For ferromagnetic nickel, it does not seem to be necessary to exploit fully these 
new fiexibilities. Therefore, we restrain ourselves to the normal-conducting state. 

2 Gutzwiller variational theory 

In this section we first define Gutzwiller variational wave functions for multi-band Hub- 
bard models. In the limit of infinite coordination number, Z — > 00, the expectation value 
of the Hamiltonian can be written in terms of an effective single-particle Hamiltonian. We 
use this expression as an approximation in our numerical studies on fee nickel {Z = 12). 

2.1 Definitions 

2.1.1 Gutzwiller-correlated wave functions 

We consider the Gutzwiller-correlated wave functions 
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Here, |\l/o) is a normalized quasi-particle vacuum so that we may later apply Wick's 
theorem. Expectation values in l^'o) are denoted by 

(i)o = (^o|i|^o) • (2) 
Expectation values in |\I^g) are denoted by 
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The Gutzwiller correlator is a product of local correlators, 

Pc^n^Glm, (4) 
m 

which induce transitions between atomic configurations on lattice site m, 

-PG|m = X! '^-fi,/2|m|A)mm(-^2| • (5) 
h,h 

Atomic configurations with at least one electron are defined by 

|/)m = n c+^a I vacuum)^ . (6) 

The operator (cm.o-) creates (annihilates) an electron on lattice site m in the local 
spin-orbitals a = 1, 2, ... , (2iVo) {No = 1, 3, 5 for s,p, d electrons). The sets / comprise all 
possible atomic occupancies with |/| = 1, 2, . . . , (2iVo) electrons. I = (0), (di), (ai, (72), . . . , 
(cTi, . . . , cr2Ar^) denotes configurations with no, one, two, . . . , 2No electrons. A fixed order 
among the orbitals is assumed, i.e., we have (Xi < cr2 < • • • < cr2No i^i the definition of the 
sets /. For example, there is only one configuration {2No) = (cti, . . . , <72No) = (!> • • • ? '^N^) 
which represents a fully occupied atom. All standard set operations are defined for the 
sets /, e.g., / = Ji U J2 is the set union, and I = Ii\ I2 denotes the set of elements in Ii 
which are not in I2. We define the complement of / as / = {2No) \ I. Altogether there 
are Dat = 2^^° atomic configurations, Dat = 4, 64, 1024 for s,p, d electrons. 

We demand that + I/2I is even in (jS)). This restriction permits pairing correla- 
tions with an even number of Fermions as in BCS theory. Furthermore, we address only 
Hermitian Gutzwiller correlators, -Pcim = -Pcim- Consequently, we demand 

Therefore, we can find atomic states |r)jn such that 

Pg\u, = Y1 -^r|m^Tir|m (8) 
r 

with real parameters Ar|m- Here, 



^r|m = |r)mm(r| 



(9) 



is the projection operator onto the state |r)m. In terms of atomic configurations we have 

|r)n. = E^/,r|m|/)m , I /)m = E ^/,r|m 1 m (10) 

/ r 
where the unitary matrix diagonahzes the matrix Am, 

Ari|m(5ri,r2 = E ^/i,ri|mA7i,/2|m^/2,r2|m ■ (H) 

In general, the state |r)ni can be chosen as the eigenstates of some effective local atomic 
Hamiltonian H^f. 

2.1.2 Multi-band Hubbard models 

In this work we focus on ferromagnetism of transition metals. Therefore, we study the 
multi-band Hubbard model for electrons on L lattice sites, 



with 



H = f + H^, (12) 



T= E ^u'ctA<.'= EeL'(k)a.^,A,.', (13) 

iJjCrcr' k;(T(T' 

C(k) = ^Ee"''^'"^^^u'' (14) 



1,J 



where we used translational invariance, tf^' = t^^' {i — j) to introduce the bare energy- 
band matrix elements e'^^,(k). Note that we set f^f = 0, i.e., the kinetic energy operator 
does not include the local potential terms ei^aa' = Co-o-' (e.g., exchange splittings, spin-orbit 
interaction) which are contained in the general, purely local Hamiltonian 

^at = E E f^/l,/2|m|/l)mm(/2| = E E ^r|m^r|m • (1^) 

m "1 r 

The projector onto the atomic eigenstate |r)in on site m is denoted mY^^. The unitary 
matrix diagonalizes the matrix 11^, 

%|m^ri,r2 = E ^h,Tr\mUhMmBj^;r2\m ■ (16) 

Il,l2 

For translation invariance, we have E^^^ = -Ef, f^/i,/2|m = f^/1,/25 ^h.hlm = ^h,i2y 
Ar|m = Ar. Note that the eigenstates \T)^ of H^t need not be identical to |r)ni. Therefore, 
the number of parameters A/^^/j in our general definition of the Gutzwiller correlator (jSJ 
is of the order of -D^^ [= O{10^) for d electrons]. In most practical applications, it is 
impossible to scan such a large variational space numerically. Therefore, one has to work 
with eigenstates |r)in of some effective atomic Hamiltonian Hlf which is characterized 
by a few physically motivated parameters. In this way, the numerical problem basically 



reduces to the minimization of the ground-state energy with respect to the real variational 
parameters Ar- This minimization problem is numerically tractable because the number 
of variational parameters is now reduced to the dimension of the atomic Hilbert space, 
Dat = C(IO^) for d electrons. 

The values for the electron-transfer parameters and the correlation parameters 
f//jj2 (or, equivalently, for e°^/(k) and 3p) depend on the specific material under investi- 
gation. We shall specify them for nickel in Section El 

2.2 Effective single-particle Hamiltonian for nickel 

For our studies on nickel we simplify the calculations by working with the correlator 

Pcim = ^ >^rrhr\m (17) 
r 

in which we choose the states |r)jn to be identical to the eigenstates of the atomic Hamil- 
tonian in (|T5|l . |r)m = |r)ni- Therefore, our main numerical problem is the minimization 
of the variational energy with respect to the parameters Ap. 

In order to determine the variational ground state of H one first needs to calculate 
the expectation value 

{H)G = E^-i{\r},{\^o)}) (18) 

of the Hamiltonian H as a. function of the variational parameters Ar and the one-particle 
wave-function I^E'o)- In the limit of infinite spatial dimensions the expectation value of 
the local Hamiltonian (fTH|) becomes 

E''\{mr}) = Lj2Ermr . (19) 
r 

The probability for finding the atomic state |r)ni, 

mr = (mr|m)G = Ap(mr|m)o , (20) 

may be used to replace the original variational parameters Ar because (mr|m)o is a simple 
function of the local density-matrix for the non-interacting system, 

Caa' = (Cm,<TCm,<7')o • (21) 

The one-particle product state |\E'o) is the ground state of the effective single-particle 
Hamiltonian 

f^^ = 5: 6..,(k)c+A,.' , (22) 

k;o'(T' 

e,,.(k) = Y.QTYAYi^)+Vaa' , (23) 

77' 

where the tensor elements Q^y are known functions of the parameters mr and the density 
matrix Co-o-', and are Lagrange parameters which act as effective on-site energies, see 



Sections ini and m The one-particle eigenvalues (= -Eicr) of T"^^ can be interpreted as quasi- 
particle energies within a Landau-Gutzwiller Fermi-liquid theory, see [HIj and Sectional 
As we will see in the case of nickel, these eigenvalues are in very good agreement with the 
band energies as found in ARPES experiments. For the kinetic energy we find 

(T(t',77' k 

(24) 

The local density matrix Co-o-', via eq. ()2H). and, therefore, the variational energy E™ 
may be considered as functions of mr and r/o-o-'- The remaining numerical task is then to 
minimize E™ with respect to these parameters. There are, however, still up to {2No)'^ + l 
constraints which have to be respect during the minimization, see (jTfj) - ()50|) . For example, 
the completeness of the states apparently leads to the constraint X^r^^r = 1- There 
is no simple recipe for the most efficient way to implement these constraints in the numer- 
ical minimization procedure. In previous works \^2^ nU\ we have proposed the following 
strategy: We perform a transformation of the local orbital basis onto a new one in which 
the local density matrix is diagonal. In this new basis, the main {2No)'^ constraints could 
then be implemented by the diagonalization of a 2A''o-dimensional matrix which we called 
the 'Z-matrix'. In the following section we present results using this strategy. The final 
numerical results should not be affected by the way the constraints are implemented. 

3 Results for ferromagnetic nickel 

In this section we first specify our model parameters. Then we give some details of 
our variational calculations. Finally, we discuss our results for the quasi-particle band 
structure and the magnetic anisotropy. 

3.1 Model specifications 

The multi-band Gutzwiller theory is not an ab-initio theory. It is based on a multi-band 
Hubbard Hamiltonian whose parameters need to be specified. For metallic systems, the 
use of a Hubbard model is well justified: screening is very efficient so that the effects 
of the Coulomb interaction between electrons at distances larger than the inverse Fermi 
wave number k^^ can be incorporated in the 'bare' band structure, and only the local, 
effective Coulomb matrix elements need to be taken into account explicitly. 

3.1.1 Electron-transfer amplitudes and local potential terms 

For nickel, we consider a minimal model which includes only those bands which are partly 
filled within a paramagnetic LDA calculation. Therefore, our multi-band Hubbard Hamil- 
tonian H comprises of 18 spin orbitals, namely 3d, 4s, and 4p. The non-magnetic local- 
density approximation to DFT provides the LDA band structure e^°^(k). We represent 
the bare energy-band matrix elements e°^/ (k) in the kinetic energy operator T in terms of 
real electron-transfer integrals t"""" (r) in the two-center approximation |^ which range up 



to third nearest neighbors, Nr < 3. We choose these parameters and the local potential 
terms eaa' in such a way that our tight-binding fit reproduces the LDA band structure. 
In this fit we include information on the symmetry of the single-particle states 
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-0.6960 


0.2300 


-0.4780 


0.3150 


-0.0481 
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-0.2092 


0.0524 


-0.0848 


0.0336 


-0.0007 
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-0.0439 


-0.0023 


-0.0245 


-0.0011 


0.0024 



Table 1: Electron-transfer parameters t„^/{r) to A^r-th neighbors (all energies in eV). 

For the optimal choice, the k-averaged root-mean-square deviation between our tight- 
binding fit and the LDA band structure up to 2eV above the Fermi energy is 15meV 
for the 3d bands and for the Asp bands. Our values of the electron-transfer parameters 
are summarized in Table ^ These model parameters have already been used for the 
calculations reported in Refs. jHOJ El] ■ 

For the on-site parameters we find for both spin species 
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7.7835 eV , 


^3d{t2g),3d{t2g) 


= ^t2g = 


-0.0290 eV , 


e3d(eg),3d(eg) 


= ^eg = 


0.0436 eV . 



where is the result from the fit to LDA and is used in practical calculations for 
a better agreement between our Gutzwiller theory and ARPES experiments for nickel, 
see below. For the total crystal-field splitting we find = et^g — = — 0.0726 eV. 

The energetically highest-lying state of pure d character is (purely t2g), 0.18 eV 
above the Fermi energy. The state X2 (purely eg) lies 0.025 eV above Ep. The width of 
the 3d bands can be estimated from -^(Xs) — E{Xi) = 4.45 eV (Xi is predominantly of 
d{eg) character), or from E^X^) — E{Li) = 4.63 eV (Li is predominantly d{t2g))- 

Next, we address the L2'-state. Experiment locates this state at about 1.0 eV below E-p. 
The SDFT calculations for ferromagnetic nickel as well as our Gutzwiller theory find 
the L2'-state about 0.3 eV below the Fermi energy when we use the LDA 4p orbital 
energy e^^"^ = 8.5335 eV. Therefore, we shift the 4p orbital energy by 0.75 eV to e^^^'^* = 
7.7835 eV, and find the L2'-state 0.97 eV below Ep. In the following we present results 
which use e^*^'^* = 7.7835 eV. This choice enhances the 4p partial density by approximately 
0.1 electron and, correspondingly, enhances the 3c? hole charge by the same amount so 
that we work with = 8.78 |21]. The remaining 1.22 valence electrons are about evenly 
split between the 4s level and the three 4p levels. 



3.1.2 Atomic interactions 

As described above, the paramagnetic LDA calculation provides a first single-particle 
contribution to the local Hamiltonian, 

-^at ^ — E ^crcrC^^cr^m,a (26) 
m,cr 

with \a) = |4s, 4j9, 3(i) ® | t)i)- addition to these single-particle contributions, we 
include the spin-orbit interaction for the 3d electrons only, \a) = \3d) | j, i). Hereby we 
restrict ourselves to the dominant, purely atomic contributions of the form 

^it^^ = E ^{(^\^x^x + ly^y+t^zW')c^,aCm,a'- (27) 

Here, ( is the strength of the spin-orbit coupling, cr^.^.z are the three Pauli matrices, and 
^x,y,z are the Cartesian components of the vector operator for the angular momentum. For 
the spin-orbit coupling constant we choose ( = 0.080 eV as in Refs. (SHI ESI- 

The second part of the local interaction is the two-particle Coulomb interaction, 

^i? = E E W(^i,^2;fX3,a4)c+_^^c+^^^c^,,3C^,<,, . (28) 

The intra-atomic interactions in the 4s and Ap shell are rather weak when compared to the 
broad Asp energy bands. Thus, we expect only small correlation effects in these bands. We 
also neglect correlations between Asp electrons and 3d electrons beyond those contained 
in the LDA band structure. This is a more serious approximation as we neglect magnetic 
polarization effects on the Asp bands and thus may underestimate the Asp contribution to 
the magnetic moment. Under these assumptions, the spin-orbit sum in ()28|) runs over the 
3d orbitals only. The Hamiltonian of the full atomic problem, Hs,t = -f^lt''^ + H^''^ + 
provides the matrix elements between two atomic configurations Ii and I2 in (jlSp . 

Ul,J, = m{h\Hil^^ + Hil'^ + ^i?V2)m . (29) 

The atomic eigenstates |r)ni can be written as product states, 

|r)m = |r3d)m|r4s)m|r4):3)m 

In our Gutzwiller theory we correlate the 3d states |r3rf)in only, i.e., the variational par- 
ameter Ar is independent of the Asp configuration |r4s)m|r4p)ni- This leads to a numeri- 
cally tractable problem with O{10^) variational parameters for the atomic 3d states. 

Naturally, we must not completely ignore the Coulomb interaction of the Asp electrons. 
In this case a big charge fiow from the 3d to the Asp bands and unphysically small 
occupations of the 3d shell would result. One way to overcome the charge fiow problem 
is to introduce a chemical potential which keeps the 3d partial charge fixed during the 
calculations, = 8.78; see jSl] for an alternate method which gives essentially the same 
results. 



(30 



Gutzwiller theory produces an optimum magnetic (spin-only) moment which is too 
large by about 10%. It should be mentioned that the SDFT results show a similarly large 
overshooting of the magnetic moment. To best compare with experiment, we thus carry 
out calculations at fixed magnetic moment for the total energy, with either the spin-only 
moment fixed to the (spin-only) experimental value of figpin = 0.55/iB [30] or with the 
total magnetic moment fixed to the full experimental moment of fi = 0.606;Ub when the 
spin-orbit coupling is included |H7j . 



3.1.3 Parameters for the Coulomb interaction 

It remains to determine the Coulomb interaction parameters W(cri, (J2; o"3, CT4) in (j2l 
In the spherical-atom approximation, there are only three independent interaction pa- 
rameters, namely the Slater-Condon integrals Fq, F2, and F4, from which all Coulomb 
interaction parameters can be determined. We prefer the Racah definition, see |3Hj, where 
the parameters A, B, and C are used which are linear combinations of the Slater-Condon 
integrals. 

The spherical-atom approximation is excellent in cubic systems. In principle, the 
Coulomb interaction among 3d{t2g) electrons may differ from the Coulomb interaction 
among 3d{eg) electrons because the radial parts of their orbital wave functions can be dif- 
ferent. Measurements of d-d transitions of magnetic impurities with cubic site symmetry 
in non-magnetic oxide hosts show that these differences are marginal. The d^ multiplets 
'^H and ^P, which are accidentally degenerate in spherical-atom approximation, split by 
cubic two-particle corrections but not by the crystal field. However, they are found to 
be degenerate within the experimental resolution of some meV jHEj- In test calculations 
we have used different A parameters for d(t2g) electrons and d{eg) electrons but we have 
not found any features in our quasi-particle energies which would indicate a failure of the 
spherical-atom approximation. 

The Racah parameters B and C are related to the Slater integrals F2 and F4. They 
determine the splitting of the multiplets of a specific 3(i" configuration. They can be 
determined experimentally from d-d transition spectra for magnetic impurity ions in non- 
magnetic insulating hosts. Typical values for all kinds of transition metal ions are tab- 
ulated in |38J- It is found experimentally that the ratio C/B varies smoothly between 
4 and 5; C/-B = 4 is obtained theoretically when hydrogen 3d wave functions are used. 
Experimentally, the values for B and C for Ni^+ and Ni^+ ions are known When 
we linearly extrapolate these values to a neutral atom we find our values for nickel as 
B = 0.09 eV and C = 0.40 eV so that we employ a ratio C/B = 4.4. The values for 
B and C are close to the 'bare' atomic values, i.e., the screening appears to be of little 
importance to B and C. For a fixed ratio C/B, we may replace the two interaction pa- 
rameters B and C by a single effective parameter J, as is quite often done in the LDA+[/ 
literature. This exchange coupling J is related to B and C by 

B C J 
2 + 5=7- 

For our nickel values for B and C we find J = 0.88 eV, a value very similar to the ones 
used by Anisimov et al. [12] and others. 



The Racah parameter A (basically the U parameter of the Hubbard model) determines 
the separation of the various multiplets. The 'bare' values, as calculated, e.g., from 
atomic wave-functions are of the order of 25 eV, as discussed already in Herring's book on 
magnetism [2^] and confirmed recently in jlHI- There is a technique to extract U param- 
eters using 'constraint' density-functional theory and supercell geometries. It has been 
found that for minimum-basis models the [/-values are smallest. For example, in the case 
of the cuprates, values ?7 < 2 eV were found for a single-band Hubbard model whereas a 
three-band model including the oxygen 2p states employs Udd ~ 8eV [H]. Similar results 
have been found in the case of BaBiOa For 3c? band models, values of [/ = 4 . . . 6 eV 
have been reported For the 4s-4p-3d multi-band case of nickel, no such calculations 
are available but values U ^ 10 eV appear to be reasonable. We choose A in such a way 
that our Gutzwiller theory reproduces the experimental 3d band width. For our 18 orbital 
basis we find A = 9eV to reproduce best the quasi-particle 3d band width W = 3.3 eV; a 
value of y4 = 12 eV leads to W = 3.0 eV. If we worked with a model of only ten 3d spin 
orbitals, the corresponding band- width reduction is achieved for A ^ 4 eV. 

3.2 Details of the calculations 

The minimization problem consists of two parts. We have to optimize the Hubbard inter- 
action through the variational parameters mr ('internal' variation) and the kinetic energy 
through the one-particle product state |\E'o) and its quasi-particle excitation spectrum 
('external' variation). During the variations we keep the partial densities fixed through a 
'chemical potential' and work for fixed magnetization (spin-only in the absence of spin- 
orbit coupling). These two minimization steps are done recursively. The separation of 
'internal' and 'external' minimization offers the advantage that computer time-consuming 
integrations in momentum space are reduced to a minimum. 

Experimentally, the anisotropy energy is -Eaniso ~ —3 /xeV per atom. Therefore, our 
total-energy calculations must reach an accuracy of (9(0.1 /leV) per atom. 

3.2.1 Band calculations 

The calculation of the kinetic energy contribution to the variational ground-state energy 
requires the sum over quasi-particle energies up to the Fermi energy Ep. In the absence of 
the spin-orbit coupling we use the point-group symmetry of the fee lattice and calculate 
all quantities in the irreducible part of the Brillouin zone. This amounts to a reduction of 
the zone volume by a factor of 48. In order to generate our mesh in momentum space, we 
divide the distance between the F and the X point by 20. This is the standard k-mesh 
used in all our calculations. For the calculations with cubic symmetry this results in 916 
points and 4.0 ■ 10'^ tetrahedra for the mapping of the Fermi surface in the irreducible 
part of the Brillouin zone. For the calculations including spin-orbit coupling, we use the 
same mesh. Now that the k-summation runs over the full Brillouin zone, we sample over 
3.4 ■ 10"^ k-points and 1.6 ■ 10^ tetrahedra. 

For the numerical integrations required for the calculation of the kinetic energy, we 
use the tetrahedron method as described in P^. Additionally, in all calculations we 



employ a refinement of tlie numerical integration near Ep wliereby we clieck wlietlier 
or not a tetraliedron contains states witliin a certain energy sliell around Ep. If so, we 
divide the tetrahedron into eight smaller ones and repeat the integration using the smaller 
tetrahedra. For the spin-only (full moment) nickel calculations our integration refinement 
near E-p results in 4.0-10^ (1.7-10®) additional tetrahedra in momentum space close to E-p. 

In order to achieve an accuracy of 10"'' eV for the ground-state energy, we need to 
suppress numerical noise. For example, a numerical integration yields the elements of 
the single-particle on-site density matrix C (12111 . For cubic symmetry, some diagonal 
elements should be identical but they show some scattering. Even worse, off-diagonal 
elements which should be zero by symmetry are finite, of the order O{10~^), albeit we 
perform our calculations in double precision. The main reason for the numerical noise 
is the limited accuracy of the tetrahedral method and of the integration refinement. In 
principle, numerical noise can be reduced by choosing a finer k-mesh; finite computational 
resources bar this route to higher accuracy. Instead, we take advantage of all symmetries 
of the problem, e.g., in a system with cubic symmetry the on-site density matrix C 
becomes diagonal. Therefore, off-diagonal matrix elements can be set to zero, and diagonal 
elements are replaced by their average. 

3.2.2 Variational procedure 

The variational energy has to be minimized with respect to the internal variational param- 
eters nip and the external parameters rj^^/. The latter can be interpreted as the effective 
on-site energies, crystal-field splittings, exchange splittings, and spin-orbit couplings. As 
shown in Section IXTl for cubic nickel, symmetry allows three external variational parame- 
ters for the 3d electrons: the effective crystal-field splitting e^f and the effective exchange 
splittings A^^(t2g) and A^^(eg), which are different for 3d(t2g) and Sd{eg) electrons. Sym- 
metry also allows exchange splittings for the Asp orbitals but these turn out to be of very 
minor importance, as long as we do not include the atomic electron-electron interaction 
between 3d electrons and Asp electrons, see Section 13.11 When we include the spin-orbit 
coupling (j27|) there is another external variational parameter, C*^^, the magnitude of the 
spin-orbit coupling for our quasi-particle description. 

The optimization of the external parameters for a given tensor Q is analogous to 
the Hartree-Fock variational theory, i.e., it is a familiar self-consistent band-structure 
problem. We would reproduce the results of a pure mean-field LDA+U approach if we 
set Q'^Y = 5^,cr'^7',cr'- In the case of nickel, however, we found that Q deviates significantly 
from this simple diagonal form. This explains why mean-field theories fail to describe the 
electronic properties of nickel correctly. 

The variation of the internal variational parameters rrip was the centerpiece of our 
programming efforts. All our schemes use random number generators for the starting 
values of the variational parameters so that we can easily test the stability of our results. 
We repeat the runs for different seeds and study the scatter of the results, not only in the 
total energy, but also in the values of the individual variational parameters. 

For fixed external variational parameters, the kinetic energy changes only due to 
the internal variational parameters mp. The tensor Q in ()24|) is a product of matrices q 



whose entries can be written in the form 



<lcTua2= S(cri,a2;ri, 12)^/^7^- (32) 

Most of the elements of the tensor E can be prepared before we start the variation of 
the internal variational parameters, i.e., they do not change during the optimization of 
the parameters mp. Note that the eigenvalues and eigenvectors of the matrix Z feed into 
the tensor S. The matrix Z of size 10 x 10 has to be recalculated and diagonalized at 
each internal variation step as it ensures fulfillment of the conditions ()47|) and ()48|). see 
Section I^?T1 However, only a few of the elements in S are affected by these recalculations. 

The present code simply varies each mr one by one, calculates the matrix Z and the 
corrections to the matrix q in each step, computes the matrices q and C, and carries out 
the total-energy summations both for the interaction energy and the kinetic energy. The 
use of the matrix q considerably speeds up the calculation and variation of Q and thus of 
the kinetic energy. 

A large number of variational parameters is very small, e.g., the probability mr for a 
local configuration Ni^"*" is tiny. We use a cut-off below which such variational parameters 
are set to zero and are no longer varied. In general, big variational parameters influence 
the energy more strongly than small ones so that we vary big variational parameters more 
often. Typically, we perform 1.0 • 10^ to 5.0 • 10^ variations to reach convergence. Hereby, 
the energy reduction in the last 2.0 ■ 10^ variations is less than 0.1 /ieV. When we repeat 
the run for different seeds, the total energy scatters by less than 1 /xeV. The individual 
scatter of band and interaction energies is of the order (9(/ieV). The scatter in the most 
significant mr values is of the order C(10~^), the same holds true for the diagonal elements 
of the matrix q. 

3.3 Quasi-particle band structure without spin-orbit coupling 

First we address the case where we set the spin-orbit coupling to zero and work with a 
fixed spin-only moment of /ispin = 0.55. Figure^and table |21 summarize our results which 
are in very good agreement with ARPES experiments. 

Our calculations show the expected dependence of the ?>d band width on the magnitude 
of the Racah parameter A via the hopping reduction tensor Q. As can be seen from tableEl 
the ?)d band width agrees very well with experiment for A = 9 eV. The Racah parameters 
only indirectly influence the quasi-particle bands through the matrix Q whose elements 
determine the reduction of the electron-transfer amplitudes. Modifications of the values 
for A, B, and C cause overall energy shifts but they only weakly influence specific band 
states. Now that the value of A is so much bigger than B and C, there is very little 
dependence of the 3d band width on the value of the Racah parameters B and C. 

A central result of our Gutzwiller theory is the anisotropic exchange splitting. The 
magnetization of our sample leads to the effective exchange splittings for 3d{eg) states 
and 3d{t2g) states, A'^^{t2g) = 0.25 eV, and A^^{eg) = 0.11 eV which differ by more than 
a factor of two. Correspondingly, in the quasi-particle band structure near Ep, states of 
pure t2g character (e.g., X5, Wi>) exhibit a large splitting of majority and minority spin 
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Figure 1: Grey-scale plot of the negative second derivative of the ARPES intensity for 
nickel with respect to energy, —(Pl/dE'^, on a logarithmic scale (insets: linear scale) for 
the TKX and TX directions of the Brillouin zone. The dispersionless structure at Ep is 
due to a residual Fermi edge resulting from indirect transitions. Some bands (A2, Xi, 
S2, S3, Si towards F) are not seen due to unfavorable matrix elements, depending on 
geometry and chosen final state |30j. Theoretical curves are G-DFT, see table |21 

states of about 400 meV whereas the states with pure Cg character (e.g., X2) exhibit small 
splittings of about 150 meV to 200 meV, see table |21 An important consequence of this 
anisotropy is the result that the state X21 lies slightly below the Fermi energy. This result 
is robust in all our calculations with Gutzwiller theory, e.g, it also applies for A = 12 eV. 

The large anisotropy of the exchange splitting has been found experimentally by spin- 
polarized ARPES studies along the A-line of the Brillouin zone jlHI, and is well reproduced 
by our Gutzwiller theory, see table El In contrast, the SDFT energy bands show a much 
larger but rather isotropic exchange splitting because the exchange-correlation potential is 
constructed both from superposed 3d(t2g) and Sd{eg) charge densities, i.e., the exchange- 
correlation potential is independent of the orbitals. In our Gutzwiller theory the two 
values for the exchange splittings result from two independent variational parameters. 

It is relatively easy to understand the reason for the large anisotropy of the exchange 
splittings. The width of the 3d{t2g) bands is mainly determined by the first-neighbor 
electron-transfer integral {dda)i, the largest 3d-3d matrix element, see tabled The width 
of the 3d{eg) bands does not depend on the {dda)i integral, but only on the much smaller 
{dda)2 integral to second nearest neighbors. In addition there is some first-neighbor 
3d{eg)-3d(t2g) hybridization via the large {dd7!')i integral. Thus, the most strongly anti- 



bonding bands are the 3d{t2g) bands which form the top of the 3d bands. Electron-electron 
interaction enforces the exchange splitting of the bands up to the point where the limit of 
strong ferromagnetism is reached and the majority spin bands are completely filled. In the 
minority spin bands, some energy can be gained when the occupation of the anti-bonding 
3d{t2g) bands is reduced, at the expense of the occupation of the top of the 3d{eg) bands, 
i.e., of the occupation of the X21 state. 
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Table 2: Binding energies in eV with respect to the Fermi energy E-p (> for occupied states). 
(. . .) indicates the spin average, error bars in the experiments without spin resolution are given 
as lb. Theoretical data show the spin average and the exchange splittings in square brackets. 
The second column denote the orbital character of the states, t = t2g, e = eg, capital letters: 
dominant character. The spin-polarized data (A^-j) were taken at fractions / of the TL distance, 
with the emphasis on the analysis of the exchange splittings. G-DFT and SDFT calculations 
are without spin-orbit coupling and a fixed spin-only moment of //spin = 0.55/iB- 



3.4 Magnetic anisotropy 



We now turn on the spin-orbit coupling. Again, we perform fixed-moment calculations, 
now with the full magnetic moment /x = 0.606 fi-g- Our total-energy calculations including 
spin-orbit coupling are two to three orders of magnitude more time consuming than the 
calculations discussed in Section l331 because our k-sums run over the full Brillouin zone, all 
hundred elements of the matrix q must be calculated and all thousand internal variational 
parameters must be optimized, and additional external variational parameters appear, 
e.g., the effective spin-orbit coupling C^^. 
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Table 3: Theoretical binding energies in eV with respect to the Fermi energy Ep (> 
for occupied states) with and without spin-orbit coupling. The second column denote 
the orbital character of the states, t = t2g, e = Cg, capital letters: dominant character. 
For the calculations without spin-orbit coupling, G-DFT, we used a spin-only moment 
of yUspin = 0.55 /iB- In the spin-orbit calculation, G-DFT(SO), we worked with a fixed 
moment of 0.606 /xb along the (lll)-direction. The directions A and the L-points in the 
Brillouin zone are not equivalent: the first and the second energy belong to the directions 
parallel and perpendicular to (111), respectively. 



Our optimal ground state displays an orbital moment of /^orb = 0.052 /jb which agrees 
very well with the experimental value of Horh = 0.0507 /^b CO- In the Hartree-Fock limit 
of our method, Ar = 1, we obtain /iorb,HF = 0.032 /^b- Obviously, the orbital-moment con- 
tribution to the total moment almost doubles when we apply our Gutzwiller theory. This 
is due to the fact that correlations induce non-diagonal elements in the matrix q; they 



become nonzero because the spin-orbit coupling lifts the cubic symmetry. An analysis 
of the quasi-particle bands shows that the non-diagonal elements of the matrix q induce 
electron transfers with spin flips in the effective single-particle Hamiltonian. These pro- 
cesses give an additional contribution to the orbital moment. The effective value for the 
spin-orbit coupling becomes smaller, C,'^^ = 0.070 eV for the optimum ground-state energy. 
This is a 12% reduction compared to the bare spin-orbit coupling constant ( = 0.080 eV. 
Calculations with a fixed = ( yield higher ground-state energies and an approximately 
10% larger orbital moment. 

The central external variational parameters, the effective exchange splittings A'^^{t2g) 
and A'^^{eg) as well as the effective crystal-field e^f, are very similar to the spin-only 
calculations. Consequently, the quasi-particle bands do not change very much from the 
quasi-particle bands as described in Section EiHl However, degeneracies at high-symmetry 
points and lines are lifted. In table El we give some of the energies at high-symmetry 
points and lines, taken from a calculation where the magnetic moment lies parallel to 
the (lll)-direction, the easy axis of nickel. Near the Fermi energy the splittings should 
be detectable by high resolution ARPES. Note that the state which derives from the 
X2;-state still lies slightly below the Fermi energy. 
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Figure 2: Cut of the Fermi surface with the plane TXK and the plane TXL in the 
Brillouin zone. Open symbols are data of Stark, communicated in [IH!, filled dots are 
de-Haas-van- Alphen data from and lines are the results from Gutzwiller theory. 

The spin-orbit coupling causes the former majority and minority spin bands to couple 
so that band crossings of spin-up bands and spin-down bands are avoided. As a conse- 
quence, the topology of the Fermi surface changes slightly. These changes affect only a 
few areas on the Fermi surface. In figure |2l we show a cut of the Fermi surface in the 
TKX plane and in the TXL plane, respectively. 

The theoretical Fermi surface agrees very well with the de-Haas-van-Alphen data 



of Stark, which have been communicated through the paper of Callaway and Wang [IH] . 
However, the experimental data exhibit oscillations in these cuts which are not reproduced 
by our Gutzwiller theory. In fact, no band calculation whatsoever has produced these 
wiggles. We suspect that these wiggles are artifacts of the reconstruction procedure from 
the raw data: in a magnetic field of several Tesla, the Fermi surface topology itself can 
change, as we now discuss. 

When we align the magnetic moment along the (OOl)-direction, the quasi-particle 
energy of the X-states at (100) and at (010) differ from the energy of the (001) X-state; 
the three states are degenerate when the moment is along the (111) direction. The spin- 
orbit coupling induces a splitting of the formerly degenerate X21 state by 30 meV, and the 
energy of the (001) X2j-state in the presence of the spin-orbit coupling lies approximately 
6 meV above the Fermi energy whereas the other two state lie 25 meV below Ep. Thereby 
our Gutzwiller theory confirms the scenario put forward by Gersdorf fF7\ . 

The differences in ground-energies of the states with the magnetic moment aligned 
in the (111) direction and with the magnetic moment aligned in the (001) direction, 
-Eaniso = -^111 — Eqqi, is -Eaniso, Gutz ~ "10 /xeV per atom. Thus we get the correct sign and 
order of magnitude of the experimental energy difference, -Eaniso = —3 /xeV per atom. In 
order to improve the confidence in our results, we have repeated our numerical integrations 
over the Brillouin zone with a denser mesh but for fixed external variational parameters. 
Our calculations using a denser mesh with up to 2.5 ■ 10^ k-points in the Brillouin zone 
confirm the result of a (111) easy axis in nickel. Therefore, we are confident that Gutzwiller 
theory calculations produce reliable results for the magnetic anisotropy energies. Further 
studies with even larger k-meshes are in progress. 



4 Diagrammatic approach 

In this section we develop a general scheme to evaluate expectation values for single-site 
and two-site operators using the variational analog of Feynman diagrams. 



4.1 Conditions 

We need to evaluate expectation values of operators Oi and Oij = cl Jc^ J for i 7^ j with 
the wave function ((T)), 
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The square of the local Gutzwiller correlator can be reduced to 

-^Glm = E '^/l,-f2|m|A)mm(/2| 



(33) 
(34) 

(35) 



with 



(36) 



We aim at a diagrammatic calculation of (Oi)G and (Oij)G- In such an approach, the 
sites m 7^ i, j in ()33|) and ()34|) will play the role of inner vertices. 

In a first step, we eliminate all local contributions at the inner vertices which arise 
after the application of Wick's theorem to the numerators in (j33|) and plj) . As in a 
previous work [22], we expand Pq\^ in terms of new local operators. 
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The Hartree-Fock operator for an even number ri > 4 of Fermion operators is defined 
recursively. 
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The prime in fl39j] indicates that 2 < 7^ < n — 2 is even. For n = 2 we have 

[aia2\^^ = (0102)0 • (41) 



By construction, the operator in (j37p generates diagrams with exactly (I/2I) non-local 
lines that enter (leave) the lattice site m. All diagrams with trivial local Hartree-Fock 
bubbles are automatically excluded. 

For our diagrammatic evaluation we demand that at least four lines meet at every 
inner vertex, 

xi,Mm = for + 1/21=2. (42) 

We can fulfill these conditions by restricting our parameters Xj-^^j^lm- Using the form ^7\i 
together with one easily sees that 
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must be fulfilled. This follows from the fact that, apart from the trivial term in (jSZ|), at 
least four lines meet at every vertex; a contraction with two 'external' fermionic opera- 
tors leaves at least one uncontracted pair whose local contraction, however, vanishes by 
construction of new local operators (jHTj) . Explicitly, 

Il,l2 
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Here, we introduced the expectation value 

rril,I,\m = (^0|/l)mm(/2|^0) 

= (^0|ri/i,/2jmri/,u/2|ml^0) , (51) 

with 

^'|m= 11(1- V<x) (52) 

and 

^m.cr = Cjjj (jCm,(T • (53) 

Equations p7j) - (jHn|) can be used to fix some of the variational parameters A/^^jim, for 
example those with |Ji| + IJ2I < 1. 

4.2 Diagrammatics 

With the help of (j37j) we expand the product over the squares of local Gutzwiller corre- 
lators in (|33|) [and in (j34|) ] in the form 



m^ij m ^ mi,m2 

where the prime on the sums indicates that all lattice sites are different from each other 
and from i [and j]. When Wick's theorem is applied to the numerators in ()33|) and in (jMj) 
we can introduce new contractions between two Fermi operators am.o- and dn,a (a = c or 
a = c^) 

In most cases, n 7^ m holds and the new definition of a contraction reduces to the 
usual one because the extra term proportional to 5n,m vanishes. In those cases where 
Fermi operators on the same site are contracted the contribution must vanish because 
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The new contraction in 



we subtracted all local Hartree contributions in 
fulfills this condition. 

The application of Wick's theorem thus gives a result which we would have obtained if 
we had worked with Grassmann operators instead of Fermion operators in the definition 
of the expectation values; all Grassmann operators af^ anti-commute with each other 



, (56) 



so that local contractions always vanish, 

(«f,A')o = 0. (57) 

The use of Grassmann operators instead of Fermi operators also shows that we may now 
drop the restrictions on the lattice sums because all contributions with two lattice sites put 
equal vanish due to the anti-commutation relation between the corresponding Grassmann 
operators. In this way we have generated a diagrammatic theory in which lines between 
two vertices n and m are given by the one-particle density matrices K^^.^^, defined 
in (jSni), and xj-^j^i^ gives the strength of a vertex (Ji,/2|m) with |Ji| in-going lines and 
I/2I out-going lines. 

Now we are in the position to apply the linked-cluster theorem so that all disconnected 
diagrams in the numerator in ()33|) [and in ()34|) ] cancel the corresponding denominator. 
Then, the calculation of the expectation values (j^Hj) and (jH^ is reduced to the sum over 
all connected diagrams according to the Feynman rules with lines and vertices as defined 
above. 



5 Exact results for infinite coordination number 

By construction, at least four lines meet at every inner vertex. Due to the absence of 
local Hartree-Fock contributions two inner vertices are always connected by at least three 
independent paths of lines (there is no diagram with a single inner vertex). For lattices 
with Z nearest neighbors, 

<.|n,.' - O ; (58) 

for a simple cubic lattice in d dimensions, Z = 2d. Therefore, the contribution of all 
diagrams with inner vertices vanishes at least proportional to 1/ \/d in the limit of infinite 
dimensions. Thus, not a single diagram needs to be calculated in infinite dimensions, i.e., 
we arrive at the result 

(Oi)G = (^ol^GiAPGiil^o) , (59) 

(6i,j)G = (v&o|PG|iPG|jOijPG|iPG|j|*o) (60) 



in (i = cxD dimensions. The result (j6(J|) can be further simplified because in ci = cxo 
dimensions there can be only a single line connecting the two external vertices i and j. 



5.1 Local quantities 

Now we address the local operators for the particle densities, 



the pairing densities, 
and the local interaction 



Hat = X! f^Jl,J2|m|-^l)mm(-^2| ■ 



In infinite dimensions we have 



Here we defined the matrix elements for the local densities as 

^hMli = X] '^^i,/2|i'^/3,-f4|ii(-^2|c[^aCi,(T'|-^3)i 

= X K,iua\-Aiua',h\i^sgn{a, /)fsgn(a', /) 

where 

fsgn((7,7) = (7Ua|c+|7) 

is ±1 if it takes an odd/even number of anti-commutations to shift the state 
proper place in the ordered sequence /. Correspondingly, we find 

Si,aa' = {Si,a(T')G = {PG\iCi,aCi,a' PG\i) 

Ii,h 

Here we defined the matrix elements for the local pairing amplitudes as 

^h,I4i = X -^^i,'f2|i'^/3,-f4|ii(-^2|Ci,aCi,o-'|-^3)i 

= X '^/i,/|i'^/u(aa'),/4|ifegn(cT, /)fsgn(cT', /)fsgn((T, a') . 
Lastly, along the same lines we find 

E^' = (^at)G 

i Ji,J2 

= E E ^Jl,-f4|i^/i,-f4|i ' 
i /l,/4 

where 

Uh,U\i = E '^/l,/2|i-^/3,/4|i^/2,/3|i • 
-f2,/3 



5.2 Single-particle density matrices 

Now we evaluate the single-particle density matrix 

= {c'l^aCi,a')G (71) 

and the one-particle pairing matrix 

Si,aja' = (Ci,crCj,o-')G (72) 

foriT^j. 

We start with Oij = c^^Cj^^^i in (jMjl and evaluate 

pG\iCt.pG\i = E C;„/,|i,J/l)ii(/4| (73) 

with 

= E -^-fi,^2|i'^/3,/4|ii(-^2|c;|'^|-^3)i 

= E fsgn(cr,/)A7,ju<7|iA/,/4|i . (74) 

In d = oo dimensions only a single line can join the two external vertices i 7^ j in (jUUI) . 
In the first case, the contraction at i is done with a creation operator c[ty. There are two 
possibihties, it can either be 7 G Ji, with 

|/i)ii(/4| = (fsgn(7, Ji \ 7) ^/A7,/4|i^/iu/4|i) ' (75) 

or 7 G Ji U I4,, with 

|/l)ii(/4| = - (fsgn(7, h) n7i,/4U7|i"'/iU/4U7|i) • (76) 

Altogether, the amplitude for such a process is given by 

^7= E Efsgn(7,^)("^5iiClu7,/'|i,a-"^?',/U7li4,/ii,,.) ■ (77) 

Here, we introduced the expectation value 

= ("'^^'|i^/U7'U7|i)o • (78) 

In previous works [2312015 the factors g?^^ were denoted by ^Jq\^. In general, however, 
g?^^ is not a positive real number. 

The second case is only possible if there are superconducting pairing correlations in 
|\E'o). Then, the line between i and j can also be an anomalous propagator, i.e., the 
contraction at i can be done with an annihilation operator Q^-y. Again, we obtain two 
contributions, 7 G /4, which gives 

|/i)ii(/4| = % (fsgn(7, h \ 7) ^/i,/4\7|i"'/iu/4|i) ' C^S) 



and 7 G /i U /4, which contributes 

|/i)ii(/4| = Ci^-y (fsgn(7, Ji) n/,U7,/4|i"'7iU74U7|i) • (80) 
Therefore, the amphtude for the anomalous process can be written as 

%ay= E Efsgn(7,^)("^?v|i4,/U7|i,<7 + "^/U7,/'|i4/'|i,-) • (^l) 
/{7^/) I' 

With these definitions we find for the one-particle density matrix 



77 



Qla'rQL'Y{cLK,y)o + gi,^7gj,a'7'(^i.7Cj,7')0 (§2) 



in the limit of infinite dimensions. 

Analogously, for the one-particle pairing matrix we find 



Si,a;}(7' — X![^i,cr7'?j,a'7'(Ci,7Cj,7')o + Q'i%7Q'j*CT'7'(Ci|7Cj^y 

77' 

in the limit of infinite dimensions. 

6 Variational ground-state energy 

In this section we specify the variational problem in the limit of high dimensions. We use 
translational invariance to render the problem numerically tractable. 

6.1 Variational problem 

The variational ground-state energy for the multi-band Hamiltonian (fT^ is given by 

^^^^({A/„/.|i},{|^o)}) = (i^)G. (84) 

With the help of the results of Section this expectation value is readily evaluated in 
infinite dimensions. We find 

^™^({A/„7.|i},{|^o)}) = i5'^^({A/,,,|i},{|^o)}) + i5^*({A,,,,|i},{|vl/o)}), (85) 
^'^°({A/„/.ii},{|^o)}) = E E C'"'^'(^o|£-X'l^o) , (86) 

i,j;(To-' a, «'=+,- 

E-\{h,j,\^,{\^o)}) = E E f^/i,/4li<,/4|i • (87) 

i hji 



where Ui-^j^\i is given in (ffn|l . The superscripts a, a' = +/— were introduced in order 
to distinguish creation and annihilation operators (c~ = c). Mathematically, these su- 
perscripts are interpreted as real numbers, e.g., c(-") is equal to c"*" or c when a = — or 
a = +, respectively. For the hopping amplitudes in ()86p we find 

rr' 

The interaction energy is given in terms of the variational parameters A/^_/2|i purely 
local expectation values in l^'o)- With the help of Wick's theorem these expectation values 
can be expressed solely in terms of the local single-particle density matrix with entries 

C^^l.' = {cl./C')o . (90) 



6.2 Minimization 



6.2.1 Effective single-particle SchrSdinger equation 

We have to optimize the variational ground-state energy (jHK|l with respect to all param- 
eters {A/ij2|i} and the set of all one-particle quasi-particle vacua {|\l/o)}, whereby the 
equations (PTjl - lj^Uj) have to be obeyed. The latter may be fulfilled by fixing the param- 
eters {A/j^/2|i} 1-^1 1 5 l-^al < 1, which are then functions of C\ and {Ai} = {\i^j2\'i} with 
l-^il; I -^2 1 > 2. We introduce the tensors r^i of Lagrange multipliers with entries ?7f^"/ and 
the multiplier E'^^ in order to fulfill equation ()90|1 and to ensure the normalization of 
In addition, we fix the average number of particles = nL in \^g) with the help of the 
Lagrange multiplier A. Then, the variational ground-state energy E^^^ is given by 



Tpvar 
^0 



^ ^Minimum ({A;}, {Q}, {r^;}. A, {|^o)}, 

{Ai},{Ci},{7yi},A,{|*o>},i?SP 



E.{...) = E^-^({\,},{c-^{\m,)})-E^''{{m,\m,)-i) 

[ctt' - (^o|a?.c?:,|vi/o)) + c.c. 



(91) 



(92) 



A(E^i,-({M,{a})-iv) . 



The minimization with respect to |\E'o) can be carried out explicitly and leads to the 
effective one-particle Schrodinger equation 



-iSP 



(93) 



I E E[te'"^' + 25i,jc:')^^X' + h-c. 

i,CT;j, cr' aa' 

^0 + ^ E E [[Cf"'^' + '^iJ - + h-c-] > (94) 



2 , , 

i,(t;j, cr' aa' 



where 



aa ' li.aa 



+ c.c. 



(95) 



We assume that the optimum wave-function |\l/o) is the ground-state of the Hamilto- 
nian (j94j) . In this way, |\l/o) becomes a function of {Ai}, {Q}, {f/i}, and the remaining 
task is to find the minimum 



TTivar 
^0 



Minimum ^^({Ai}, {Q}, {^i}, A) , 

{Ai},{Q},{r,i},A 

({Ai},{a},{r/i}) + i?^*({Ai},{Q}) 



E 



(96) 
(97) 



i.aa Lace 



Y: K:.'C^,:.' + c-c. - S.a'A K.. {Ai}, {Q} 



n 



6.2.2 Translational invariance 

To make further progress we assume translational invariance. Then, the single-particle 
Hamiltonian in the Bloch basis with wave vectors k has the form 



(Q:k),(TC(-a'k),cr' ) 



(98) 



kicrcr' aa' 



where the coefficients 



1.J 



Vera' 



'(-o.'){-a) 
'la' a 



'laa' 



(99) 



are the elements of the matrices eg"'. In (I99|) we used the bare energy-band matrix (jl4p 
and introduced the coefficients 



Qtr' 



— \aa 
tt' 



* — * — — 

Ira'ir'a' ~ Qra^r'a' ' Vtt' 



— l-'^aa' * — * — * * ^ ; 

Vrr' ~ Iralr'a' ^ Qralr'a' 5 Vtt' 



Qtt' 
Vrr' 



(100) 
(101) 



Let us regard the operators eg ^ as elements of vectors eg. Then, the effective Hamiltonian 
T^^ can be written in terms of matrix products. 



The matrix in (jl02p may be diagonalized by means of a Bogoliubov transformation 



where we introduced new fermionic operators and matrices u'^' with elements h'^^ and 
u""'(k), respectively. For reasons of consistency the matrices v!^' obey the symmetries 

utt = {u^-Y , «+k = (%^)* • (104) 

The fermionic commutation rules of the new operators are ensured when the matrix 
in p03|) is unitary. In ()102|) - ()104|) we have used the standard notations for the transpo- 
sition M^^ , the complex conjugate M*, and the Hermitian conjugate of a matrix M. 
After the diagonalization T^^ becomes 

k 

= To + Eo + Y,E^.fi+^h^,^ , (105) 

k 

where the real quantities i?k,7 = S^yE^ylh) are the elements of the diagonal matrix 

^k = 2E(4""^")^C'Sk", (106) 

aa' 

and 

^o = -^ETr(^k) . (107) 
^ k 

For the derivation of ()105|) - ()106|) we have used the following symmetries of the matrices 

^k ; 

= -(?^'")^ = -(?^'")* («^a'), (108) 
e-'k" = - (i^'Y = - {a^a') (109) 

which follow from (jU^ and 

4,7(k) = «7'(k))* = <7'(-k) (110) 

which result from hermiticity of T and the fact that our electron-transfer amplitudes are 
real. Note that the matrices £t^'" , e^" , the operators h^,^, and the energies i?k,7 still 
depend on the parameters A, rj, C. The Fermi-gas ground state of ()1()5|1 is given as 

l^op,^,^)) = n'/^k,7|vacuum) . (Ill) 

k,7 

Here, the prime indicates that only those single-particle states with 

^k,7(A,^,^) <^F = (112) 

are filled and 

(a, v, C) = To (r}) + Eo (A, v,C) E^,, (A, v, C) (113) 

k,7 



is the energy of the 'Fermi-gas' ground state. In principle, the variational ground-state 
energy can now be calculated by a minimization of the energy-functional 



era' act' 

+LA[Y,n^J\,C)-n] (114) 



with respect to all parameters Xj in A, C, rj, A, i.e., 

Eo™'- = E,(AO,C°,r7°,A°) (115) 

where 

^E, = . (116) 

Minimizing E^. in this straightforward way requires rather time-consuming calculations of 
expectation values in |\l/o) for every single variation of our parameters. Such a strategy is 
prohibitive, in particular due to the large number of variational parameters A. Therefore, 
one needs to develop more sophisticated numerical strategies when our general theory is 
to be applied to reahstic model system. 



7 Landau— Gutzwiller quasi-particles 

The Gutzwiller theory, as described in the previous section, is a variational approach and 
by itself a method to examine ground-state properties only. As shown in Refs. ji31"| ITT] , it 
is possible to use the approximate Gutzwiller-correlated ground state for the calculation 
of elementary excitations. In particular, the 'band structure' E^^^ derived in the previous 
section gives the variational excitation energies of Landau-Gutzwiller quasi-particles. 



7.1 Definition of single-particle excitations 

Originally, Gutzwiller-correlated wave functions were introduced as Fermi-liquid ground 
states for ordinary metals jS]. More generally, they are Fermi-liquid ground states in the 
sense of an adiabatic continuity to some non-interacting reference system which can be a 
Fermi gas, a BCS superconductor, or a band insulator. For all such cases, quasi-particle 
and quasi-hole states are readily identified. To this end, we define quasi-particle and 
quasi-hole creation operators as 

:= PgV(^g)-^ 
These operators obey usual Fermi anti-commutation relations, 

[K,r^^p'y]+ = Sp,p' Stt' , (119) 
[K,r^K',r']+ = [Vp,r,Vp',M = 0, (120) 



(117) 
(118) 



and they create quasi-particles/quasi-holes in the variational Gutzwiller ground-state, as 
can be seen from 

e;^,v^,,\^G) = e(EF-i?p,.)|^G), (121) 

Vptrl^c) = QiE^,r-EF)\^G). (122) 
We define the quasi-particle or quasi-hole excitation energy as 

E±(p, r) = ± (i^r (P, - E^,^^) - /x± (iVr (P, r)-N) , (123) 

where the upper and lower sign corresponds to quasi-particles and quasi-holes, respec- 
tively. Here, 

are the expectation values of the energy and the average particle number in the quasi- 
particle/quasi-hole states 

I^Si:^) = e+JvPG) = PG/i^,.|^o)=PG|^iT^), (125) 
1^^'^) = V^vl/^) = PcVrl^o) = Pcl^i^V (126) 

The one-particle state |\l/o) is defined according to (jlllj) in terms of the operators 
Note that |\E'o) and /ij^,^ have to be used with their optimum values, given by the pa- 
rameters A°, r/'', C°, A°. The quasi-particle excitation energy ()123p is measured from the 
(variational) chemical potential /i^ which describes the energy for adding a particle to 
the system or subtracting it, respectively, 

/i± = ±{El^\N ± 1) - ^r(^)) • (127) 

Eq'^'^{N) is the (variational) ground-state energy for a system with N electrons. For a 
metallic or superconducting system, /i+ = /i~ = /i. For an insulating system. A™ = 
^'^ — fi~ defines the gap in our variational theory. The variational ground-state energy 
flllSp depends on the particle density n = N/ L implicitly, due to the n-dependence of all 
parameters in X^,C^,rf,A^, and, explicitly, through the term —An in ()114|) . Therefore, 
/i± is given by 



L dn ^ L dxi 



J-A»^-A». (128) 



where we used (lll(jp and the fact that, in an insulating system, both ^± and A*^ are 
discontinuous as a function of n. 



7.2 Quasi-particle dispersion 

We write the energies -E™'' and -E±'^(p, 7) in p2H|l as 
1 



771 var 



(T CT'', 77', QQ;' 



nat/TO /=;0 



where 



(129) 



(130) 



The tensor A^'^ gives the change of C due to the creation of a quasi-particle or a quasi- 
hole, 

A^^^,(p, r) = ± (y,7'(p, r) - <"(p, r)) , (131) 

where we introduced 

v;7'(P'^) = «"^-(p))*<-(p). (132) 

The sum over k in (jl29p is readily evaluated, 
E^L'(k)(^S^^|fi(„k).7C^Kk),yl^o^^) = KZ'''''C^',C',v') (133) 

k 

± (k,7'(p, r)e^(p) - 17;"(p, r)6°,^(p)) . 
To order (1/L)° the expansion of 6E^^^ = 2{E^'{p, r) - E^''') with respect to A^'^ yields 



^^at 



+ E 

pp' ,aa 

+ E 



crcr 

era' a/^aa' 



C=CO pp 



OQ 

-cr',77' I'-'pp' 



+ 



C=CO 



A^^^,(p,r) (134) 
(A^^p'(P,r); 



C=CO. 



We use 



±^p,. = ±^ E E«'""'fe''(P'^)^^(p)-^7';"(P'^)^^(p) 



2 ^ ^ 

(TO-', 77' aa' 



+^ E [((C'T-(<'7)V±'Sv(P'^)+c-^- 

77',aQ!' 



(135) 



to evaluate the first line in ()134p . For the second and the third line we employ 



'■77' dQTJ'''' 



d 



c=c° 



-T 



'Off 



^^PP' 



(136) 



C=CO 



and the corresponding equation where C"^' is replaced by (Cp^"') • In we have used 

the fact that the dependence of K^^/ '^^ on C can be ignored in the derivative of E^^ 
with respect to C because l^'o) is an exact eigenstate of T"^^ and E^^ is its eigenvalue. 
Using ()114|1 and ()116|) we can rewrite this expression as 



dE 



SP 



pp 



dE 



at 



C=CO 



C=CO 



+K(«T-("?';)°)-vAii: 



(137) 



C=CO 



where, again, the same equations hold with C and t] replaced by C* and r^*, respectively. 
With these results and 



ivr(p,r)-iv=5: 



A 



±;pp' 



, + C.C. 



(138) 



we find i?-|-(p, 7) = E^^^ for the quasi-particle and quasi- hole excitation energy. This 
relations shows that the 'band structure' E-p^.^ describes meaningful quasi-particle bands. 



8 Outlook 

Our Gutzwiller theory provides a very good description of the quasi-particle band struc- 
ture of ferromagnetic nickel. We obtain the correct exchange splittings and, in particular, 
we reproduce the experimental Fermi-surface topology. We find the correct (lll)-direction 
of the magnetic easy axis and the right order of magnitude of the magnetic anisotropy. Our 
theory even reproduces the experimentally observed change of the Fermi-surface topology 
when the magnetic moment is oriented along the (001) axis. 

Our investigations show that SDFT should be improved along the following lines. 
First, it needs to incorporate an orbital dependence of the exchange-correlation potential 
in order to reproduce the anisotropy of the exchange splittings. It may also be used 
to correct the partial densities of the "id and 4s]9 electrons. Second, it should allow an 
effective spin-orbit coupling. 

Thus far, our calculations have focused on the quasi-particle band structure at zero 
temperature. Several extensions of the Gutzwiller theory are possible. First, as shown 
in 1:4 7j, we can also address magnetic excitations. Therefore, Gutzwiller theory should 
reproduce the spin-wave spectrum of nickel. Second, we can generalize the idea of ele- 
mentary excitations to transitions between atomic states which are induced by photoe- 
mission. In this way we should also be able to make contact with the experimentally 
observed '6eV-peak' [iH] . 

We emphasize that our Gutzwiller theory as presented here is not limited to nickel. 
We may apply it equally to other transition metals and their compounds, for example 
NiO jini- Our theory also covers superconducting pairing which we may apply to super- 
conductivity in multi-band systems. The Gutzwiller approach can and will be applied to 
many interesting correlated-electron problems in the future. 
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